In [88]:
import numpy as np
%matplotlib inline
import matplotlib.pyplot as plt
from scipy.integrate import odeint
from scipy import signal
Modèle sous forme d'espace d'état :
$\begin{cases} \frac{d}{dt}\begin{bmatrix}i(t)\\ \omega_{R}(t) \end{bmatrix} & =\begin{bmatrix}-\frac{R}{L} & -\frac{K}{L}\\ \frac{K}{J} & -\frac{f}{J} \end{bmatrix}\begin{bmatrix}i(t)\\ \omega_{R}(t) \end{bmatrix}+\begin{bmatrix}\frac{1}{L} & 0\\ 0 & -\frac{1}{J} \end{bmatrix}\begin{bmatrix}V_{S}(t)\\ \tau_{L}(t) \end{bmatrix}\\ \begin{bmatrix}i_{obs}(t)\\ \omega_{Robs}(t)\\ V_{Lobs}(t) \end{bmatrix} & =\begin{bmatrix}1 & 0\\ 0 & 1\\ -R & -K \end{bmatrix}\begin{bmatrix}i(t)\\ \omega_{R}(t) \end{bmatrix}+\begin{bmatrix}0 & 0\\ 0 & 0\\ 1 & 0 \end{bmatrix}\begin{bmatrix}V_{S}(t)\\ \tau_{L}(t) \end{bmatrix} \end{cases}$
In [89]:
def heaviside(x):
if x == 0:
return 0.5
return 0 if x < 0 else 1
In [110]:
@interact
def _(x0 = slider(0,10,10/100.0),
x1 = slider(0,1,1/100.0),
sim_duration = slider(10,30,default=10,label="Simulation duration"),
sim_nb_samples = slider(10,100, default=100,label="Samples number"),
steptime = slider(0,10,default=5,label="Step time"),
stepgain = slider(-5,5,default=2,label="Step gain"),
R = slider(0,2,2/100.0,default=1),
K = slider(0,0.02,0.02/100.0,default=0.01),
L = slider(0.5/100,1,1/100.0,default=0.5),
J = slider(0.01/100.0,0.02,0.02/100.0,default=0.01),
f = slider(0,0.2,0.2/100.0,default=0.1),
action=selector(["MCC outputs", "RRA\'s"],label="Simulation choice")):
##### Modèle oracle
A = matrix([[ -R/L, -K/L ],
[ K/J , -f/J ]])
B = matrix([[ 1/L , 0 ],
[ 0 , -1/J ]])
C = matrix([[ 1 , 0 ],
[ 0 , 1 ],
[ -R, -K ]])
D = matrix([[ 0, 0 ],
[ 0, 0 ],
[ 1, 0 ]])
##### Modèle défectueux
Af = matrix([[ -R*1.3/L, -K/L+1 ],
[ K/J , -f/J ]])
Bf = matrix([[ 1.2/L , 0 ],
[ 0 , -1/J ]])
Cf = matrix([[ 1 , 0 ],
[ 0 , 1 ],
[ -R, -K ]])
Df = matrix([[ 0, 0 ],
[ 0, 0 ],
[ 1, 0 ]])
ss_mcc = signal.lti(A,B,C,D)
ss_faulty_mcc = signal.lti(Af,Bf,Cf,Df)
t = np.linspace(0,sim_duration,sim_nb_samples)
u1 = [stepgain*heaviside(x-steptime) for x in np.nditer(t)]
u2 = np.array(np.zeros_like(t))
U=np.vstack((u1,u2)).transpose()
X0 = (x0,x1)
tout, Y, X = signal.lsim(ss_mcc, U, t, X0)
tout, Yobs, Xobs = signal.lsim(ss_faulty_mcc, U, t, X0)
# Nice examples of interactive calculus : https://wiki.sagemath.org/interact/calculus
# Plot help : http://sage.brandoncurtis.com/plotting.html
# Plot help2 : http://doc.sagemath.org/html/en/reference/plotting/sage/plot/plot.html
# Matplotlib examples : http://matplotlib.org/examples/index.html
if action == "MCC outputs":
f, axarr = plt.subplots(3, sharex=True, figsize=(10,6))
axarr[0].plot(t, Y[:,0], 'b')
axarr[0].grid(True)
axarr[0].set_ylabel('$i_{obs}(t)$',fontsize=20)
axarr[1].plot(t, Y[:,1], 'r')
axarr[1].grid(True)
axarr[1].set_ylabel('$\omega_{R_{obs}}(t)$',fontsize=20)
axarr[2].plot(t, Y[:,2], 'g')
axarr[2].grid(True)
axarr[2].set_ylabel('$V_{L_{obs}}(t)$',fontsize=20)
plt.xlabel('Time (s)',fontsize=16)
f.subplots_adjust(hspace=0.3)
plt.show()
elif action == "RRA\'s":
n = len(t)-2
U1 = np.zeros((n, 2))
U2 = np.zeros((n, 2))
Yobs1 = np.zeros((n, 3))
Yobs2 = np.zeros((n, 3))
Yobs3 = np.zeros((n, 3))
r_ar1 = np.zeros(n)
r_ar2 = np.zeros(n)
r_ar3 = np.zeros(n)
r_ir = np.zeros(n)
r_1 = np.zeros_like(t)
r_2 = np.zeros_like(t)
r_3 = np.zeros_like(t)
##### Calcul des RRA's auto-redondantes et inter-redondante
for i in range(len(t)-2):
U1[i] = U[i,:]
U2[i] = U[i+1,:]
Yobs1[i] = Y[i,:]
Yobs2[i] = Y[i+1,:]
Yobs3[i] = Y[i+2,:]
##### Calcul des valeurs des RRA's
r_ar1[i] = J*L*Yobs3[i,0] - J*U2[i,0] - U1[i,0]*f + (K^2 + R*f)*Yobs1[i,0] + (J*R + L*f)*Yobs2[i,0] - K*U1[i,1]
r_ar2[i] = J*L*Yobs3[i,1] - K*U1[i,0] + (K^2 + R*f)*Yobs1[i,1] + (J*R + L*f)*Yobs2[i,1] + R*U1[i,1] + L*U2[i,1]
r_ar3[i] = K^2*Yobs1[i,2] + J*R*Yobs2[i,2] + J*L*Yobs3[i,2] + K^2*U1[i,0] + J*R*U2[i,0] - K*L*U2[i,1] + (R*Yobs1[i,2] + L*Yobs2[i,2] + R*U1[i,0])*f
r_ir[i] = J*L*Yobs2[i,2] + J*R*U1[i,0] - K*L*U1[i,1] + (J*R + L*f)*Yobs1[i,1] + (K^2*L + L*R*f)*Yobs1[i,0]
##### Calcul des RRA's directes
for i in range(len(t)):
r_1[i] = R*Yobs[i,0] + K*Yobs[i,1] + Yobs[i,2] - U[i,0]
r_2[i] = Yobs[i,0] - Y[i,0]
r_3[i] = Yobs[i,1] - Y[i,1]
f, axarr = plt.subplots(7, sharex=True, figsize=(10,12))
axarr[0].plot(t[0:-2], r_ar1, 'b')
axarr[0].grid(True)
axarr[0].set_ylabel('$r_{ar1}(t)$',fontsize=20)
axarr[1].plot(t[0:-2], r_ar2, 'r')
axarr[1].grid(True)
axarr[1].set_ylabel('$r_{ar2}(t)$',fontsize=20)
axarr[2].plot(t[0:-2], r_ar3, 'g')
axarr[2].grid(True)
axarr[2].set_ylabel('$r_{ar3}(t)$',fontsize=20)
axarr[3].plot(t[0:-2], r_ir, 'g')
axarr[3].grid(True)
axarr[3].set_ylabel('$r_{ir}(t)$',fontsize=20)
axarr[4].plot(t, r_1, 'b')
axarr[4].grid(True)
axarr[4].set_ylabel('$r_{1}(t)$',fontsize=20)
axarr[5].plot(t, r_2, 'r')
axarr[5].grid(True)
axarr[5].set_ylabel('$r_{2}(t)$',fontsize=20)
axarr[6].plot(t, r_3, 'g')
axarr[6].grid(True)
axarr[6].set_ylabel('$r_{3}(t)$',fontsize=20)
plt.xlabel('Time (s)',fontsize=16)
f.subplots_adjust(hspace=0.3)
plt.show()
In [ ]: